Day 5 – Spatial data lab

SICSS 2025

During the lab, you will need to load some libraries. You can install and load them with the following code:

list.of.packages <- c("tidyverse", "osmdata", "sf", "tmap", "leaflet")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

#Data manipulation
library(tidyverse)

#Spatial data
library(osmdata)
library(sf)
library(tmap)
library(leaflet)

In this lab, we are going to see how to handle spatial data using OpenStreetMap.

OpenStreetMap (OSM) is a free, open geographic database updated and maintained by a community of volunteers via open collaboration 1.

It provides an easily accessible API from which we can retrieve spatial data. In R, we can use the osmdata package to get tidy output from the API query.

Retrieving data on OSM

With osmdata, a query can be done like this:

  1. Define an area from which you want to collect spatial information. This area is called a bounding box. Use opq and/or getbb.

  2. Add the features you want to collect. A list can be found here (or use osmdata::available_features()). A feature consists of a key (a broad category, e.g. a building) and a value (e.g., a church). Use add_osm_feature.

  3. Convert the result to a sf object. Use osmdata_sf.

For instance, let’s find out the schools of Norrköping:

OSM data can be a bit messy, since this is a collaborative effort. Here, schools are sometimes points, lines, polygons, or multipolygons.

Mapping data

It is a bit uncertain as to what we captured with this query: are the schools really in Norrköping? Are they really schools? It’s always better to double-check with a map!

Here, there are two options: static or interactive maps. Several packages allow to map sf objects. We are going to focus on ggplot2 for static maps and tmap for interactive ones. leaflet is also a good options for interactive maps.

library(tmap)
tmap_mode("view") #Set interactive mode

#Create a data.frame with all schools. Careful, this data.frame includes multiple geometries!
nkpg_schools_sf <- dplyr::bind_rows(
  st_make_valid(nkpg_schools$osm_multipolygons),
  nkpg_schools$osm_polygons[!is.na(nkpg_schools$osm_polygons$name), ],
  nkpg_schools$osm_points[!is.na(nkpg_schools$osm_points$name), ]
  ) 

#Map
tm_basemap("Esri.WorldGrayCanvas") +
  tm_shape(nkpg_schools_sf) +
  tm_polygons(popup.vars = "name") +
  tm_shape(nkpg_schools_sf[st_dimension(nkpg_schools_sf) == 0, ]) +
  tm_dots(popup.vars = "name")

Exercise 1

  1. Using the pipeline described above, find all the hospitals in Stockholm.
  2. Map the hospitals.
sklm_hospitals <- 
  #Create the query and set the bounding box
  opq(bbox = "Stockholm") %>% 
  #Select the features
  add_osm_feature(key = "amenity", value = "hospital") %>%
  #Convert to osmdata_sf
  osmdata_sf()

sklm_hospitals_sf <- bind_rows(
  st_make_valid(sklm_hospitals$osm_multipolygons),
  sklm_hospitals$osm_polygons[!is.na(sklm_hospitals$osm_polygons$name), ],
  sklm_hospitals$osm_points[!is.na(sklm_hospitals$osm_points$name), ]
  ) 

tm_basemap("Esri.WorldGrayCanvas") +
  tm_shape(sklm_hospitals_sf) +
  tm_polygons(popup.vars = "name") +
  tm_shape(sklm_hospitals_sf[st_dimension(sklm_hospitals_sf) == 0, ]) +
  tm_dots(popup.vars = "name")

Exercise 2

Using osmdata::getbb, get the boundaries of the city of Paris, and of Sweden. Check ?getbb.

#Paris 
paris <- getbb("Paris, France", format_out = 'sf_polygon', featuretype = "settlement")
plot(st_geometry(paris))

#Sweden
sweden <- getbb("Sweden", format_out = 'sf_polygon', featuretype = "country")
plot(sweden$multipolygon)